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Abstract 

The paper is devoted to the numerical solution of elastoplastic constitutive initial value prob¬ 
lems. An improved form of the implicit return-mapping scheme for nonsmooth yield surfaces is 
proposed that systematically builds on a subdifferential formulation of the flow rule. The main 
advantage of this approach is that the treatment of singular points, such as apices or edges at 
which the flow direction is multivalued involves only a uniquely defined set of non-linear equa¬ 
tions, similarly to smooth yield surfaces. This paper (PART I) is focused on isotropic models 
containing: a) yield surfaces with one or two apices (singular points) laying on the hydrostatic 
axis; b ) plastic pseudo-potentials that are independent of the Lode angle; c) nonlinear isotropic 
hardening (optionally). It is shown that for some models the improved integration scheme also 
enables to a priori decide about a type of the return and investigate existence, uniqueness and 
semismoothness of discretized constitutive operators in implicit form. Further, the semismooth 
Newton method is introduced to solve incremental boundary-value problems. The paper also 
contains numerical examples related to slope stability with available Matlab implementation. 

Keywords: clastoplasticity, nonsmooth yield surface, multivalued flow direction, implicit return- 
mapping scheme, semismooth Newton method, limit analysis 

1 Introduction 

The paper is devoted to the numerical solution of small-strain quasi-static elastoplastic problems. 
Such a problem consists of the constitutive initial value problem (CIVP) and the balance equation rep¬ 
resenting the principle of virtual work. A broadly exploited and universal numerical/computational 
concept includes the following steps: 

(а) time-discretization of CIVP leading to an incremental constitutive problem; 

(б) derivation of the constitutive and consistent tangent operators; 

(c) substitution of the constitutive (stress-strain) operator into the balance equation leading to the 
incremental boundary value problem in terms of displacements; 

(d) finite element discretization and derivation of a system of nonlinear equations; 

(e) solving the system using a nonsmooth variant of the Newton method. 
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CIVP satisfies thermodynamical laws and usually involves internal variables such as plastic strains 
or hardening parameters. Several integration schemes for numerical solution of CIVP were suggested. 
For their overview, we refer, e.g., [9] ED ED US] and references introduced therein. If the implicit or 
trapezoidal Euler method is used then the incremental constitutive problem is solved by the elastic 
predictor/plastic corrector method. The plastic correction leads to the return-mapping scheme. We 
distinguish, e.g., implicit, trapezoidal or midpoint return-mappings depending on a chosen time- 
discretization m Chapter 7]. 

In this paper, we assume that the plastic flow direction is generated by the plastic potential, g. 
If g is smooth then the corresponding plastic flow direction is uniquely determined by the derivative 
of g and consequently, the plastic flow rule reads as follows, e.g. HD Chapter 8]: 


kP = \ dg{tr,A) 

da 


g = g(a, A). 


( 1 . 1 ) 


Here, e p , A, a , and A denotes the plastic strain rate, the plastic multiplier rate, the stress tensor 
and the hardening thermodynamical forces, respectively. The corresponding return-mapping scheme 
is relatively straightforward and leads to solving a system of nonlinear equations. A difficulty arises 
when g is nonsmooth. Mostly, it happens if the yield surface contains singular points, such as apices 
or edges. Then the function g is rather pseudo-potential than potential and its derivative need not 
exist everywhere. In such a case, the rule (1.1) is usually completed by some additive formulas 
depending on particular cases of g and a in an ad-hoc manner. For example, the implementation of 
the Mohr-Coulomb model reported in [EL; Chapter 6, 8] employs one, two, or six plastic multipliers 
A, depending on the location of a on the yield surface. Since the stress tensor a is unknown in CIVP 
one must blindly guess its right location. Moreover, for each tested location, one must usually solve 
an auxilliary system of nonlinear equations whose solvability is not guaranteed in general. These 
facts are evident drawbacks of the current return-mapping schemes. 

In associative plasticity, it is well-known that the plastic flow rule (1.1) together with a hardening 
law and loading/unloading conditions can be equivalently replaced by the principle of maximum 
plastic dissipation within the constitutive model. This alternative formulation of CIVP does not 
require special treatment for nonsmooth g and enables to solve CIVP by techniques based on math¬ 
ematical programming In particular, if the implicit or trapezoidal Euler method is used 

then the incremental constitutive problem can be interpreted by a certain kind of the closest-point 
projection PEDES]. For some nonassociative models, CIVP can be re-formulated using a theory 
of bipotentials that leads to new numerical schemes Pill HE ED]. These alternative definitions of 
the flow rule enable a variational re-formulation of the initial boundary value elastoplastic problem. 
Consequently, solvability of this problem can be investigated (see, e.g., [T8, EU). Therefore, the 
corresponding numerical techniques are usually also correct from the mathematical point of view. 
On the other hand, such a numerical treatment is not so universal and its implementation is more 
involved/too complex in comparison with standard procedures of computational inelasticity. 

The approach pursued in this paper builds on the sub different al formulation of the plastic flow 
rule, e.g. HD Section 6.3.9], 

e p G \d a g(a, A) (1.2) 

for nonsmooth g. Here, d a g(a, A) denotes the subdifferential of g at (cr, A) with respect to the stress 
variable. If g is convex at least in vicinity of the yield surface then this definition is justified, e.g., by 
[ 30; Corollary 23.7.1] and is valid even when g is not smooth at a. On the first sight, it seems that 
(1.2) is not convenient for numerical treatment due to the presence of the multivalued flow direction. 
The main goal of this paper is to show that the opposite is true, by demonstrating that the implicit 
return-mapping scheme based on (1.2) leads to solving a just one system of nonlinear equations 
regardless whether the unknown stress tensor lies on the smooth portion of the yield surface or not 
at least for a wide class of models with nonsmooth plastic pseudo-potentials. Using this technique, we 
eliminate the blind guessing and thus considerably simplify the solution scheme. Moreover, the new 
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technique enables to investigate some useful properties of the constitutive operator, like uniqueness 
or semismoothness, that are not obvious for the current technique. 


1.1 Basic idea 

First of all, we illustrate the new technique on a simple 2D projective problem that mimics the 
structure of an incremental elastoplastic constitutive problem. Consider the convex set 

B := {w = (w 1 ,w 2 ) G M 2 | f(w) < 0}, f(w) := w 1 + \w 2 \ - 1, 

and define the projection w* G B of a point 2 = {z\,z 2 ) G R 2 as follows: 

IIz - w*|| 2 = min ||2 - m|| 2 , || 2 || 2 := 2 2 + 2 2 

W£B 1 

The scheme of the projection is depicted in Figure [lj 



Figure 1: Scheme of the projection. 

Clearly, the function / is convex in R 2 , nondifferentiable at w = (tfi,0) and 


v/M = i, 


w 2 

\w 2 \ 


Vm = (wi,w 2 ), w 2 ^ 0. 


(1.3) 


If 2 G B then w* = 11^(2) = 2 . Conversely, if 2 ^ B it follows from the Karush-Kuhn-Tucker 
conditions and (1.3) that the projective problem can be written as follows: find w* = 


*\T 


and the Lagrange multiplier A > 0: 

Z\ — w\ = A, Z 2 — W 2 G \d\wl\, w\ + | lUj | — 1 = 0, 

where 

Q| *i _ / i W 2 /1^2 I}i W 2 ~f~ 0) 

1 2l “ 1 [- 1 , 1 ], «$ = 0 , 


(1.4) 


To find a solution to (1.4), it is crucial to rewrite the inclusion (1.4) 2 as an equation. Observe that 

z 2 — G \d\w 2 \ if and only if w’ = (|^ 2 | — A) + —^r, 

\ Z 2\ 

where (-) + denotes a positive part of a function. This crucial transformation will be derived in detail 
in Section [ 3 ] on an analogous elastoplastic example. Thus (1.4) leads to the following system of 
equations: 

(1.5) 


l-A)+ ri’ 

* 1 1 * 1 

W 1 + | w 2 \ 

\Z2\ 
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Since (1.5) 2 implies \w^\ = (\z 2 \ — A) + , the system of three nonlinear equations reduces to a single 
one 

Z\ — A + (| ^ 2 1 — A) + — 1 = 0. 

Consequently, A can be found in the closed form as 


^ —Zi — 1 +-(—z i + \ z 2 \ + l) — <j 2 Zl - i 


+ | Z 2 | — 1 ), Z\ 


C2 


Zl - \z 2 \ 


1 < 0 , 
1 > 0 


from which one can easily compute w * = {w^w^f by (|l.5[)i and (|l.5|)2- 


1.2 Content of the the paper 

The presented idea is systematically extended on some elastoplastic models. This paper, PART I, 
is focused on isotropic models containing: a) yield surfaces with one or two apices (singular points) 
laying on the hydrostatic axis; b ) plastic pseudo-potentials that are independent of the Lode angle; 
c) nonlinear isotropic hardening (optionally). Such models are usually formulated by the Haigh- 
Westergaard coordinates. Further, the implicit Euler discretization of CIVP is considered and thus 
two types of return on the yield surface within the plastic correction are distinguished: ( i ) return to 
the smooth portion of the yield surface; (ii) return to the apex (apices). 

The paper is organized as follows. Section[2]contains some preliminaries related to invariants of the 
stress tensor and semismooth functions. Section [3] is devoted to the Drucker-Prager model including 
the nonlinear isotropic hardening. Although the plastic corrector cannot be found in closed form, 
the new technique enables to a priori decide about the return type and prove existence, uniqueness 
and semismoothness of the implicit constitutive operator. The consistent tangent operator is also 
introduced. In Section [4| we derive similar results for the perfect plastic part of the Jirasek-Grassl 
model |jT5j . In Section [5j the new technique is extended on an abstract model written by the Haigh- 
Westergaard coordinates. I 11 particular, within the plastic correction, we formulate a unique system 
of nonlinear equations which is common for the both type of the return. It can lead to a more correct 
and/or simpler solution scheme in comparison with the current technique. Section [6] is devoted to 
numerical realization of the incremental boundary value elastoplastic problem using the semismooth 
Newton method. I 11 Section [7j illustrative numerical examples related to a slope stability benchmark 
are considered. Here, limit load is analyzed by an incremental method depending on a mesh type 
and mesh density for the Drucker-Prager and Jirasek-Grassl models. 

Within this paper, second order tensors, matrices and vectors are denoted by bold letters. As 
usual, small letters are used for vectors and capitals for matrices (see Section [6]). Further, the fourth 
order tensors are denoted by capital blackboard letters, e.g., O e or I dev . The symbol ® means the 
tensor product muni- We also use the following notation: M + := {z G M; z > 0} and for the 
space of symmetric, second order tensors. 


2 Preliminaries 

2.1 Invariants of a stress tensor and their derivatives 

Consider a stress tensor a e M;/// and its splitting into the volumetric and deviatoric parts: 

& — pi + s, p: = p(cr) = ^/ : cr, s = s(<t ) = l dev : cr = cr - ^(/ : a )/. 

Here, I, ldev, V an d s denote the identity second order tensor, the fourth order deviatoric projection 
tensor, the hydrostatic pressure, and the deviatoric stress, respectively. The Haigh-Westergaard 
coordinates are created by the invariants p, g and 9 , where 

Q ■= e{v) = Ys : s = HI, 
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9 := 6(a) = - arccos 
o 


' 3 ^ J 3 


J. 


3/2 


Ms) = 


s = -s : s = -g 
2 2 


Ms) = 3 s 3 


s s cr 


Clearly, g > 0 and 9 G [0, 7t/ 3]. Since 6 is not well-defined when g = 0, the Lode angle is included in 
elastoplastic models indirectly. Usually, it is considered another invariant in the form 


g := g(cr) = gr(cos6), 


( 2 . 1 ) 


where f(-) is assumed to be a smooth function such that g (-) is at least continuous. However, as we 
will see below, it is more convenient to assume strong semismoothness of p(-). As a particular case, 
it will be considered the invariant 


Qe '■= QeW) = gr e (cos9), 


( 2 . 2 ) 


where 


r e (cos 6) 


4(1 — e 2 ) cos 2 9 + (2e — l) 2 
2(1 — e 2 ) cos 9 + (2e — 1) y/4(l — e 2 ) cos 2 9 + 5e 2 — 4e 


(2.3) 


The function r e was proposed in [35] and contains the excentricity parameter e G [0.5,1]. It holds: 
a) r e (cos0(.)) is a bounded and smooth function for any cr G M 3 *^, g(a) > 0; b ) g e (cr) = 0 when 
g = 0; c) g e = g, r e (cos9) = 1 when e = 1. 

We will also use the following derivatives: 


dp 


I ds 
3 ’ dcr 


I dev ) ^ (hr) 


dg s 
dcr g ’ 


— = - (l dev - n <g) n ), 

dcr g 


(2.4) 


86 

for 


Vg 

psin(36 ) ) 


[(n (8) n 3 )/ 


Idei;(^ 2 )] , 


dr e 

da 


—r' e (cos9) sin 9 


d9 

8a 


(2.5) 


Notice that the derivatives of g, n, 9 and r e do not exist when g = 0. Further, 9 is not differentiable 
when cr satisfies either 9 = 0 or 9 = n/3. On the other hand, r e has derivatives for such stresses [39j . 
For purposes of this paper, it is crucial to derive the subdifferential of g at cr when g(cr) = 0: 


dg(cr) = {n G Mj„ 3 | g(r) > g(cr) + n : (r - c r) Vr G M 3 ^} 

= {he M 3 ym I \\s(r)\\ > (h : I)(p(r) - p(c r)) + h : s(r) Vr G M 3 ^} 

= {he Mj„ 3 | I :h = 0, ||s(r)|| > h : s(r) Vr G M 3 ^} 

= {he Mj„ 3 | I :h = 0 , ||n|| < 1 } if g(cr) = 0 . 


If g(cr) > 0 then dg(cr) 


(n(cr)} by (2.4) 3 . It is 


readily seen that 


I : h = 0 Vn G dg(cr), 


( 2 . 6 ) 


(2.7) 


regardless g(cr) = 0 or not. 

2.2 Semismooth functions 

Semismoothness was originally introduced by Mifflin [25j for functionals. Qi and J. Sun |28] extended 
the definition of semismoothness to vector-valued functions to investigate the superlinear convergence 
of the Newton method. We introduce a definition of strongly semismooth functions [34,i EDA 1251 • 4o 
this end, consider finite dimensional spaces A" and Y with the norms ||.||x and ||.||y, respectively. 
In the context of this paper, the abstract spaces A", Y represent either subspaces of M n or the space 
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Definition 2.1. Let F : X —>• Y be locally Lipschitz function in a neighborhood of some x E X and 
dF(x) denote the generalized Jacobian in the Clarke sense |2J/. We say that F is strongly semismooth 
at x if 

(i) F is directionally differentiable at x, 

(ii) for any h G X, h —> 0, and any V G dF(x + h), 

F(x + h)~ F(x) - VIi = 0(\\h\\ 2 x ). (2.8) 


Notice that the estimate (2.8) is called the quadratic approximate property in |34] or the strong 


G-semismoothness in [HI Eg. In literature, there exists several equivalent definitions of strongly 
semismooth functions, see 


For example the condition (2.8) can be replaced with 
F'(x + h;h)- F'(x; h) = 0{\\h\\ 2 x ) V/i G X, h -G 0, x + h G V F , 


(2.9) 


where Dp is the subset of X, F is Frechet differentiable and F'(x\ h ) denotes the directional derivative 
of F at a point x and a direction h. We say that F : X —> Y is strongly semismooth on an open set 
O C A" if F is strongly semismooth at every point of O. 


Since it is difficult to straightforwardly prove (2.8) or (2.9), we summarize several auxilliary 
results. Firstly, piecewise smooth (PC 1 ) functions with Lipschitz continuous derivatives of se¬ 
lected functions belong among strongly semismooth functions H2t Especially, we mention the 
max-function, £ : M. —y R, £(s) = max{0, s} = s + . Further, scalar product, sum, compositions 
of strongly semismooth functions are strongly semismooth. Finally, we will use the following version 
of the implicit function theorem [HJ [Ml . 25] . 

Theorem 2.1. LetX : Y x X X be a locally Lipschitz function in a neighborhood of (y,x), which 
solve X(y,x) = 0. Let d x X(y,x) denote the generalized derivatives of X at (y,x) with respect to the 
variables x. If d x X(y,x) is of maximal rank, i.e. the following implication holds, 


X° x /\x = t), X°Ed x X(y,x) 


Ax = 0, 


( 2 . 10 ) 


then there exists an open neighborhood Oy of y and a function F : Oy —>■ X such that F is locally 
Lipschitz continuous in Oy, F(y) = x and for every y in Oy, X (y,F(y)) = 0. 

Moreover, ifX is strongly semismooth at (y,x), then F strongly semismooth at y. 

The semismoothness of constitutive operators in elastoplasticity has been studied e.g. in [0j, 16j 
I5TL 35l [5Ti . [6j. Namely in [Ml 36], one can find an abstract framework how to investigate it for 
operators in an implicit form. However, the framework cannot be straightforwardly used for models 
investigated in this paper. Therefore, we introduce the following auxilliary result. 

Proposition 2.1. Let p(-), s(-), g(-), d(-), n(-), r e (-) and g e (-) be the functions introduced in Section 


any cr G g(c r) = 0. Define, 

g(a)r e (cos9(cr)), g(c r) ^ 0, . , __ f p(tr)I + ^(cr)n(cr), g(c r) ^ 0, 

0, g(cr) = 0, ’ ' \ p(cr)I , q((t) = 0. 

Then the functions g, g e , g e and S are strongly semismooth in 



2.1 


Further, ler p, g 


d3x3 

"syrn 




be strongly semismooth functions and assume that g vanishes for 
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Proof. Since the functions n(-) and r e (cos #(•)) are bounded and have Lipschitz continuous derivatives 
in {er G | g(cr) ^ 0}, it is easy to see that the functions g, g e , g e and S are locally Lipschitz 

continuous in and strongly semismooth for any er G K'j) , g[a) ^ 0. Therefore, it remains to 
show strong semismoothness at er G g(cr) = 0. To this end, we show that (2.9) holds for g, g e , 

g e and S at such cr. Let r G be such that g(r) ^ 0. Then 

g(cr + er) = £q(t), n(er + er) = n(r), 9(cr + er) = 9(t) Ve > 0. 

Hence g'(cr + r; r) — g'(cr ; r) = 0, £/ e (cr + r; r) — ^(cr; r) = 0 and 
p / e (o- + r;r)-p / e (<T;r) = [p'(cr + r; r) - g'(o~, r)]r e (cos 0(t))=O(\\t || 2 ), 

S'(a + r; r) - S'(cr; r) = [p'(cr + r; r) - p(cr; r)]J + [^(<r + r; r) - ^'(cr; r)]n(r)=0(||r || 2 ), 


since the functions p, g satisfy (2.9) by the assumption. 


□ 


The function S introduced in Proposition 2.1 has the same scheme as a mapping between trial 
and unknown stress tensors for models introduced in Section [3]-[5| Here, the trial stress is represented 
by cr and the unknown stress is in the form cr = p(cr)I + g(cr)n(cr). Therefore, it is sufficient to 
prove only semismoothness of the scalar functions p, g representing invariants of the unknown stress 


tensor. The semismoothness of g e has been derived to prove Theorem 4.3 


3 The Drucker-Prager model 

3.1 Constitutive initial value problem 

We consider the elastoplastic problem containing the Drucker-Prager criterion, a nonassociative 
plastic flow rule and a nonlinear isotropic hardening. Within a thermodynamical framework with 
internal variables, we introduce the corresponding constitutive initial value problem, see m ■■ 

1 . Additive decomposition of the infinitesimal strain tensor e on elastic and plastic parts: 

£ = £ e + £ P . (3.1) 

2. Linear isotropic elastic law between the stress and the elastic strain: 

cr = D e : £ e = K(I : £ e )I + 2GI dev : £ e , O e = KI ® I + 2GI dev , (3.2) 

where K , G > 0 denotes the bulk, and shear moduli, respectively. 

3. Non-linear isotropic hardening: 

k = H{eP). (3.3) 

Here e 1 * G M + denotes an isotropic (scalar) hardening variable, n G M + is the corresponding 
thermodynamical force and H : R + —>■ M + is a nondecreasing, strongly semismooth function 
satisfying H( 0) = 0. 

4. Drucker-Prager yield function: 


f(<r,K) = /(p( cr ),^( cr ),«) 


g + gp — f(c 0 + k). 


(3.4) 


Here, the parameters r/,f>0 are usually calculated from the friction angle using a sufficient 
approximation of the Mohr-Coulomb yield surface and Cq > 0 denotes the initial cohesion. 
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5. Plastic pseudo-potential. 


30) = g(p 0)> eO)) = Y 2 1 - + vp - 

Here fj > 0 denotes a parameter depending on the dilatancy angle. 
6. Nonassociative plastic flow rule 


(3.5) 


£ p G Xdg(cr), (3.6) 

where A > 0 is a multiplier and dg(cr) denotes the sub differential of the convex function g at 


cr. Using (2.4), (2.6) and (3.5), the flow rule (3.6) can be written as 


® p = A I V I 1 I > nede(tr) 


Consequently by (3.1), (3.2) and fh7), 

& = D e : (e — e p ) = D e : e — A (' G\[2h + Kfjl 

7. Associative hardening law: 


i° = -pJflvd = k , 

OK 


8. Loading/unloding criterion: 


(3.7) 


(3,8) 


(3.9) 


5 > 0, /(<r, ft)<0, A/(<t,k) = 0. 


(3.10) 


Then the elastoplastic constitutive initial value problem reads as follows: Given the history of 
the strain tensor e = s(t), t G [to,t max ], and the initial values £ p (t 0 ) = e p , ^(fo) = e p . Find 
(<j(t)iE p {t),£ p (t)) such that 

a = O e : (e - e p ), 

& = D e : e — A (GV2h + Kfjl ) , n G dg(cr), 

e p = k, 

A>0, /(cr, H{£?)) < 0, A/(cr — 0. 



hold for each instant t G /oAmax]- 


3.2 Implicit Euler discretization of Cl VP 

We discretize CIVP using the implicit Euler method. To this end we assume a partition 

0 to ^ t\ . . . <C t k <C . . . <C tjv t max 


of the pseudo-time interval and fix a step h. For the sake of brevity, we omit the index k and write 
cr \= cr(tk), £ e(tfe), £ p £ p {t k ) and e 9 = Further, we define the following trial variables: 

£ p,tr : = £P(tk-i), £ e,tr = s(tk) — £ p (tk- 1 ) and cr tr := D e : £ e ’ tr . Then the discrete elastoplastic 
constitutive problem for the A;-step reads as follows: Given £, cr tr and e p,tr . Find cr. e? and A A 
satisfying: 

cr = cr tr — AA (G\/ 2 h + Kfjl ) , n G dg(cr), 1 

eP — £P> tr + AA£, > (3.12) 

AA > 0, /(<x, H(e*>)) < 0, AA f(tr, H(e*)) = 0. J 


Notice that the remaining input parameter for the next step, 
formula £ p = £ — 


K 1 


\t k ), can be computed using the 


: cr after finding a solution to problem (3.12). 













3.3 Solution of the incremental problem 


We standardly use the elastic predictor/plastic corrector method for solving (3.12). 

Elastic predictor applies when 

f{(T tr ,H(£ P ' tr )) < 0. 

Then the triplet 


(T = (T 


tr 


= g pfT , AA = 0 


(3.13) 

(3.14) 


is the solution to (3.12). 

Plastic corrector applies when ( 3.13| ) does not hold. Then AA > 0 and (3.12) reduces into 

AA (G\/2ri +Kfjl) , hedg(cr), 
eP = ePtr + AA£, 
f{<r,H{e?)) = 0. 


CT = <7 


tr 


(3.15) 


Since the functions / and g depend on a only through the variables g and p, it is natural to reduce a 
number of uknowns in problem (3.15). To this end, we split (3.15)i into the deviatoric and volumetric 
parts: 


s = s tr — AXGV^n, h e dg(cr), 
p = p tr — A A Kfj, 


(3.16) 

(3.17) 


where s , p tr denotes the deviatoric stress, and the hydrostatic stress related to cr tr , respectively. 
Using (2.4) 3 , the equality (3.16) yields 


s tr = 


1 , aags/A 

e ) 


s if q > 0 , 
AAG\/ 2 n if 0 = 0. 


(3.18) 


Denote g tr = ||s <r || and recall that ||n|| < 1 for g = 0 by (2.6). Then from (3.18) we obtain 

g = g tr — AAGa /2 if g > 0, 

0 > g tr - AAGa /2 if 0 = 0. 


Following the arguments developed in Section 1.1, we now rewrite (3.19) as follows: 

0 = (g tr - A A GV2) + = max ( 0 ; g tr - AAGv^) . 


Notice that (3.19) and (3.20) are equivalent. Further from (3.18)i, we standardly have: 

tv 
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Atr 


S S tr ’ r a 

n ----- — n if 0 > 0 . 

e e 


(3.19) 


(3.20) 


(3.21) 


The following theorem summarizes and completes the derived results. 


Theorem 3.1. Let f(cr tr ) H (£ p ' tr )) > 0. If (cr,e p , AA) is a solution to problem (3.15) and p = I : 
<r/3, s = I dev : cr, 0 = || s|| 7 then (p, g, e p 1 AA) is a solution to the following system: 

p = p tr — AA Kg, 

0 = (g tr - AAGV2) + , . , , 

gp = £Ptr + AA ^ I 

Sip, g,H(er)) = 0. 


Conversely, if(p,g,e p , AA) is a solution to (3.22) then {er,£ p , AA) is the solution to (3.15) where 

y tr - AA {GV2n tr + Kfjl) if g tr > AAGa/ 2 , 

(p tr - AXKfj) I if g tr < AXGy/2. 


(T = 


(3.23) 
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Notice that the knowledge of the subdifferential of g enables us to formulate the plastic corrector 
problem as a unique system of nonlinear equations in comparison to the current technique introduced 
in m Moreover, one can eliminate the unknowns p, g, eP similarly as for the current return-mapping 


scheme of this model. Inserting of (3.22) x _ 3 into (3.22 ) 4 leads to the nonlinear equation q tr ( AA) = 0 
where 

fr(7) - q(r,p‘ r , e‘ r , <? p ’ ,r ) = \J\ [f - igN) + + - 7 l<n) - £ (c„ + + 7 £)), 7 e R+, 

(3.24) 



> 0 then AA > g tr /Gy/ 2 and g = 0. 


Proof. From (3.24) and the assumptions on H, it is readily seen that qt r is a continuous and decreasing 


function. Further, q tr ( 7 ) —> —00 as 7 —* +cx) and q tr ( 0) = / (cr tr , H[e p,tr )) > 0. Therefore, the 
equation q tr { AA) = 0 has just one solution in M + . If q tr (g tr /Gy/‘. 2) < 0 then AA G (0, g tr /Gy/2). 
Otherwise, AA > g tr /Gy/2. The rest of the proof follows from Theorem 


3.1 


and the elastic prediction. 

□ 


The second part of Theorem 3A is very useful from the computational point of view: one can 
a priori decide whether return to the smooth portion of the yield surface happens or not. This is 
the main difference in comparison with the current return-mapping scheme. The improved return¬ 
mapping scheme reads as follows. 

Return to the smooth portion 

1. Necessary and sufficient condition : q tr ( 0) > 0 and q tr (g tr /Gy/2) < 0. 

2 . Find AA e (0, g tr /Gy/2): 

^ ( g tr - AA Gy/2) + 77 (p tr - AA Kfj) -f(c 0 + H(^ tr + AAO) = 0. (3.25) 


3. Set 


a = a tr - A A (Gy/2n tr + Kpl) , 


£P = £P’ tr + AAe 


(3.26) 


Return to the apex 

1. Necessary and sufficient condition : q tr [g tr /Gy/ 2) > 0. 

2 . Find AA > g tr /Gy/2 : 

rj (p tr - AA Kfj) -f(c 0 + H(e p ’ tr + AAO) = 0 . (3.27) 


3. Set 


= (p tr - AA Kfj) /, 


s p = e p ’ tr + AAf. 


(3.28) 


Nonlinear equations (3.25) and ( 3.27[ ) can be solve d by the N ewton method. Then it is natural to 
use the initial choice AA° = 0, A A 0 = g tr /Gy/2 for (3.25), and (3.27), respectiv ely. In case of pe rfect 
plasticity, H = 0, or linear hardening, Fh{£ p ) = He?, H = const., equations (3.25) and (3.27) are 
linear and thus A A can be found in the closed form. 
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3.4 Stress-strain and consistent tangent operators 


Solving the problem (3.12), we obtain a nonlinear and implicit operator between the stress tensor, 


a = a(tk), and the strain tensor, e = e{t k ). The stress-strain operator, T, also depends on e p {t k _f) 
and through the trial variables. To emphasize this fact we write a := T(e; e p (t k ~i), ^(tk-i))- 


From the results introduced in Section 3.3, we have T(e\e p (t k _i),£ p (t k _i)) = S (a tr , e p,tr ), where 


-tr 


S(a t \e p ' tr ) = 


r tr 


if q tr ( 0) < 0, 

A A ( G\f2n tr + Kfjl ) if q tr ( 0) > 0, q tr (g tr /Gy/2) < 0, 
C P tr ~ AA Kfj) I if q tr (g^/Gyfi) > 0, 


(3.29) 


where AA is the solution to (3.25), ( |3.27| ) in ( |3.29[ )2, and ( |3.29[ )3, respectively, i.e. q tr ( AA) = 0. 


Theorem 3.3. The function T is strongly semismooth in 


d3x3 
^sym 


with respect to e. 


Proof. We use the framework introduced in Section 2.2 Consider the function AA := AX(a tr , e p ’ tr ) 


satisfying AA = 0 if q tr ( 0) < 0, otherwise q tr ( AA) = 0. Applying Theorem 2.1 on the implicit function 
q tr , one can easily hnd that the function AA is strongly semismooth. Consequently, the functions 
p(a tr ,e p,tr ) — p tr — (AX) + Kfj and g{a tr , e p,tr ) = ( p tr — (AA) + Gv / 2) + are strongly semismooth. Since 
S(cr tr ,e p,tr ) = p(a tr ,e p,tr )I + g(cr tr , e p,tr )n(cr tr ), we obtain strong semismoothness of the functions 
S and T using Proposition 2.1| □ 


Notice that T is not smooth if q tr { 0) = 0 or q tr (^ r /Gy / 2) = 0 or if H has not derivative at 
e p,tr + AA£. We introduce the derivative dcr/de under the assumption that any of these conditions 
does not hold. Set H\ := H'(£ p,tr + AA£). Using (2.4), (2.5), ( |T2 ) and the chain rule, we obtain the 
following auxiliary derivatives: 


da 


tr 


dp 


tr 


= Ki¬ 


ds 


tr 


= 2G\ 


dev-) 


dg 


tr 


de 


<~v ■ o 

de de de 

We distinguish three possible cases: 

1. Let q tr ( 0) < 0 (elastic response). Then clearly, 

ckr 

de 


= 2 Gn 


tr 


dn tr _ 2 G f 
~de~ ~ ~gA “ 


n 


tr 


n tr ). 


(3.30) 


2. Let q tr { 0) > 0 and q tr (g tr /G\/ 2) < 0 (return to the smooth surface ). Then the derivative of 
(13.261) reads 


^ = De - AA^^ (I dev - n tr ® n tr ) - (G-/2n tr + Kql) ® 
de g tr de 

Applying the implicit function theorem on ( |3.25[ ), we obtain 

dAX _ G\[2n tr + qKI 
de G + Krjfj + tf 2 Hi 


Hence, 


da_ 

de 


D e - AA 


2G 2 y/2 

g tr 


(I dev - n tr ® n tr ) - {GV2n tr + Kfjl) <g) 


G\f2n tr + qKI 

G + Kqq + eHi 


(3.31) 
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3. Let q tr ( K Q tr /G\J 2) > 0 (return to the apex). Then the derivative of (3.28) yields 




I — rjl 


<9AA\ 

~dT J 


Applying the implicit function theorem on (3.27), we obtain 

<9 A A 


r/K 


de 


Hence, 




Krjfj + i 2 Hi 
Krjf) 

Krfrj + i 


I. 


I® I. 


(3.32) 


The derivatives (3.30)—(3.32) define the consistent tangent operator, T°. It is readily seen that the 


tangent operator is symmetric if ?/ = rj, i.e. for the associative plasticity. For purposes of Section [6j 
it is useful to extend the definition of T° for nondifferential points. For example, one can write 


if 


qtr{ 0 ) < 0 , 


T°(£;e p (4-i),A , (4- 1 )) = { (3A1) if q tr ( 0) > 0, q tr (g tr /G\/2) < 0, 

(3732) if qtr (e*/Gy/2) > 0, 


(3.33) 


where Hi in (3.31), (3.32) is the derivative from left of H at e p,tr + AA£. Notice that 


T°(e;e p (4_i),£ p (4_i)) e <9 e T(e; e p (4-i),e p (4-i))- 


4 A simplified version of the Jirasek-Grassl model 

The Jirasek-Grassl model was introduced in pT5j. It is a plastic-damage model proposed for complex 
modelling of concrete failure. The model has been further developed. For example, Unteregger and 
Hofstetter [38] have improved a hardening law and used the model in rock mechanics. For the sake 
of simplicity, we only consider a perfect plastic part of this model to illustrate the suggested idea and 
improve the implicit return-mapping scheme. The whole plastic part of the Jirasek-Grassl model can 
be included to an abstract model studied in the next section. 


4.1 Constitutive problem and its solution 

The perfect plastic model contains the yield function proposed in 

/(» = /(p( cr ), ^( cr )> 6*{tr)) = ^ + m o 


V6. fc fc ] ’ 


(4.1) 


where mo is the friction parameter and f c is the uniaxial compressive strength. The invariants p, g 
and g e = gr e (cos9) were introduced in Section|2j Notice that the couple ( p a , g a ) = (f c /rrio, 0) defines 
the apex of the yields surface generated by the function /. The yield surface is not smooth only at 
this apex. Scheme of the yield surface can be found in |15| . 

Further, the following plastic pseudo-potential is considered 


i ay a+y> 


(4.2) 


where 


_ P~ft /3 

rrig(p) = A g B g f c e , 


An, B g , fc, ft > 0 . 
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(4.3) 



















The subdifferential of g consists of the following directions: 


-g v {p(cr), g{cr))I + g e {p{cr), g(cr))h , h E dg(cr), 


where dg(cr) is defined by ( 2 . 6 ) and 

dl rn' g (p ) 

dp 


dg 3g mo 


P-It/3 


gviP , e) ■= 47 = 9 e (P> 0) '= = if + ~7=T, m W = A g e Ss/c • 


fc 


dg ff a/6 f c 


The fc-step of the incremental constitutive problem received by the implicit Euler method reads 
as follows. Given e := e(tk) and cr tr := O e : (s(t k ) — e v (t k _i)). Find a = cr(t k ) and AA satisfying: 


<7 = CT 


tr 


AA 


K 


SM/ 4 . 2G ( — + m ° | n 
/c J + l A + aA/ c 1 n 


, n e 9^(cr), 


AA > 0 , f(p(cr),g(cr),ge(cr)) < 0 , A\f(p(cr), g(a), g e {cr)) = 0 . 


(4.4) 


We solve this problem again by the elastic predictor/plastic corrector method. Within the plastic 
correction, we define the trial variables s tr , p tr , £> tr , n ir = s tr /g tr , 9 tr , rf = r e (cos9 tr ) and n/ 
associated with er tr and obtain the following result. 


Theorem 4.1. Let f(p tr , g tr , £>/) > 0. If (er, AA) is a solution to problem (4-4) an d p — I : erf 3, 
s = : <r ; ^ = ||s|| ; i/ien (p, £>, AA) is a solution to the following system: 


g = 


p = pP r -A\K^, 

Jc 

e ,r - A « G (I + a 


f(p,Q,er f) = o. 


(4.5) 


Conversely, if (p, g, AA) is a solution to (4-5) then (er, AA) is i/ie solution to ( 4-4) where 


CT = 


p/Apii 1 ’’ !/ ^ > AA2G (| + , 

pj */ ^ < AA2G (f + j/. 


(4.6) 


Proof. To prove Theorem 4.1 we use the same technique as in Section 3.3 It is based on the splitting 
the stress tensor on the deviatoric and volumetric parts, and on using linear dependence between s 
and s tr to reduce a number of unknowns. In particular, we have 


s tr = [ 1 + AA 2 G [ At + 


m o 


fc \/6 fj g 


for g > 0. Consequently, we obtain ( 4 . 5 ) 2 , n = n tr and also 9 = 9 tr for g > 0 using ( 2 . 1 ). Finally, 
notice that grf —> 0 as g tr —» 0 + . Indeed, £> —> 0 as —> 0+ and the function r e (cos(.)) is 
bounded. □ 

Analogously to the Drucker-Prager model, one can analyze existence and uniqueness of a solution 


to problem (4.5), and a priori decide whether the return to the smooth portion of the yield surface 


happens or not. To this end, we define implicit functions p tr '■ 7 p 7 and g tr : 7 1-4 such that 

- 7 2g(% + AtY 


+ 7 /I 




fc 


-p tr = 0 , A- 


tr 


fc VQfcJ 


= 0 , 


respectively, for any 7 > 0. The following lemma is a consequence of the implicit function theorem. 
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Lemma 4.1. The functions p tr and g tr are well-defined in M + . Further, p tr is smooth and decreasing 


in M +7 Qtr is decreasing in the interval 


0 , 


Ve fce tr 

2Gtyiq 


and its closed form reads as follows: 


Qtr{ 7) = 


i + 7f 


tv 

Q -7 


2Grri( 

VQfc 


V 7 > 0. 


(4.7) 


Now, consider the function q tr ( 7 ) := qfi/]p tr , Q tr ^ 


Qtr( 7) = f {Ptr (7)) ftr(7)» Md)^) = ^ 


+ m 0 


etr{l)rf , PM 


Vbfc 


+ 


fc 


- 1 . 


(4.8) 


Theorem 4.2. Let f(p tr , g tr , of ) > 0. Then there exists a unique solution, AA > 0, of the equation 
q tr { AA) = 0. Furthermore, problems J f.5 ) and (4-4) have also unique solutions. 

In addition, if qtr ( \/f>f c Q tr / 2Gmf) < 0 i/ien AA G (0, \/P>f c g tr /2Gmfi) and g > 0. Conversely, if 
qtr (VP>fcQ tr /2Gm 0 ) > 0 i/ien AA > \/6f c g tr /2Gm 0 and g = 0 . 


Proof. Since g tr > 0 and g' tr < 0 in 


0 , j, the functions g tr , (fi r are decreasing in this interval. 

it is follows that 


2Gmo ’ ^ iese functions vanish. Therefore, from (4.8) and Lemma 


4.1 


For 7 > 

q tr is a continuous and decreasing function in M + . Furthermore, q tr ( 7 ) —> — 00 " as 7 —» +00 and 
g(0) = f(p tr , g tr , of) > 0. Hence, the equation q tr ( AA) = 0 has a unique solution in M + . If 
qtr {VP>f c g tr /2.Gmfij < 0 then AA G ( 0 , V6f c g tr /2Gm 0 ). Otherwise, AA > V6.f c g tr /2Grrio. The rest 
of the proof follows from Theorem 4.1 and the elastic prediction. □ 


Although the function q tr is implicit the decision criterion introduced in Theorem 4.2 can be 
found in closed form. 


Lemma 4.2. 


qtr 


VQfc 


tr 


2 Gmn 


>0 if and only if p‘ 


tr 


V6K rn'Jp 0 


_tr 


2 G 


m 0 


p a > 0 , 


p a = 


= —. (4.9) 

m 0 


Proof. Since g tr ) = °> 


Qtr 


V 6 f c g t r 
2Gmn 


g \ ffill mo 


fc 


-Ptr 


VQfcf 

2 Gmn 


- 1 > 0 . 


Hence, p crit := pt 


crit _~ ( vT/eg" 


r 1 2 Gm. 0 


> p a . Using the definitions of p tr and m g , we have 


0 = p‘ 


crit 


VQg 


tr 


, V 60 1 


tr 


2Gm~“" 9 


Km'(p crit ) - p tr > p a - -—fi—Km'(p a ) - p tr . 

2Gmt] 9 


□ 


By Theorem 4.1 and Theorem 4.2 the return-mapping scheme reads as follows. 

Return to the smooth portion 

1. Necessary and sufficient condition: q tr { 0) > 0 and q tr (y/6f c g tr /2Gm 0 ) < 0. 
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2. Find p E M, g > 0 and A A E ( 0 , V6f c g tr /2Gm 0 )\ 


3. 


p + AXK 1 ^ - p* r = 0 , 


£ +AA2G(|§ + 


raq_ 

x/ftfc 


rdr 


1 ( T ) + m O 


gr l 


_L P. _ I 

Vefc ^ fc 1 


0 , > 

0 . 


(4.10) 


cr — pi F gn tr . 


(4.11) 


Return to the apex 

1. Necessary and sufficient condition : q tr (y/6f c g tr /2Gm 0 ) > 0. 

2 . Set 

P = ~, Q = 0 , <r=pl, A\ = -^—(p tr -p). 
m 0 Km' g [p) 


(4.12) 


The system (4.10) of nonlinear equations can be solved by the Newton method with the initial 
choice p° = p tr , = g tr , AA° = 0. It was shown that the system has a unique solution subject to 
< 7 tr( 0 ) > 0 and q tr (x /6 f c g tr / 2Gmo) < 0. Without these conditions, one cannot guarantee existence 
and uniqueness of the solution to (|4. 10). 


4.2 Stress-strain and consistent tangent operators 


Solving the problem (4.4), we obtain a nonlinear and implicit operator between the stress tensor, 


<r = cr(tfc), and the strain tensor, e = e(£*,). The stress-strain operator, T, also depends on £ p {tk- 1 ) 
through the trial stress. To emphasize this fact we write <j := T(e; e p {t] t _i)). We have 


T(e;eP(t fc _i)) = 



qtr{0) < 0, 

g'tr(O) > 0, q tr (V6f c g tr /2Gm 0 ) < 0, 
qtr {VdfcQ tr /2Gm 0 ) > 0, 


(4.13) 


where the function q tr is defined by (4.8) and p, g are components of the solution to (4.10). 


Theorem 4.3. The function T is strongly semismooth in 


d3x3 
^sym 


with respect to e. 


Proof. Consider the function AA = AA(er tr ) satisfying AA = 0 if q tr ( 0) < 0, otherwise q tr ( AA) = 0. 
To apply Theorem 2.1 on the implicit function q tr { 7 ) := < 7 ( 7 ; p tr , g tr ), it is necessary to show that q 


is strongly semismooth w.r.t. the variables q,p , £r r . This follows from (4.7) and Proposition 2.1 


The rest of the proof coincides with the proof of Theorem 3.3 


□l 


If q tr ( 0) = 0 or q tr ( \/6f c g tr /2Gm 0 ) = 0 then T is not smooth. We derive the derivative dcr/de 
under the assumption that any of these conditions does not hold. If q tr { 0) < 0 (clastic response) 
then dcr/de = O e . If q tr [y/6f c g tr /2Gm 0 ) > 0 (return to the apex) then dcr/de = O (vanishes). 

Let q tr ( 0) > 0 and q tr (\fQf c g tr /2Gm<f) < 0, i.e., return to the smooth portion happens. Then 
the derivative dcr/de can be found as follows. 


1 . Find the solution (p, g, AA) to (4.10). 
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2. Use (2.4), (2.5), (3.2) and the chain rule and compute: 


ds ds OS 


dO 


tr 


ds Q 

„tr 


= (hev - n tr ® n tr ), 


ds g tr sin(3# tr ) 
3. Compute: 


/ # \ / 1 + AXK 7 <M 

Jc 


c ?r y \/f\ r)r tr r)f) tr 

S [(n tr < 8 ) ( n tr ) 3 )I — I de „(ri tr ) 2 ] , = —r' e (cosO tr ) sin#* 7 ' 


<9e 


ds 


as 

de 

as 

dA\ 

as 


OZAA 

\ as / 


K 


V 


0 

mo 

fc 


1 + AA§£ 2 G(^ + 


m 'g(p) 

fc 


\ - 1 / <9p rr 


3g . mpr " 
fc V&fc 


3 e I m 0 
j? sAfc 

0 


7 


V 


,ir 

as~ 

dg tr 

as 

mpg dr 1 ? 
V6fc d£ 


(4.14) 


Notice that the matrix in (4.14) arises from linearization of (4.10) around the solution (p, g, AA). 
The matrix is invertible since its determinant is negative. 


4. Compute 


dcr dp dg dn tr 

— = I®/ + n ir <8>-^ + . 

ds ds ds ds 


(4.15) 


For numerical purposes, we use the following generalized consistent tangent operator: 

O e if q tr ( 0 ) < 0 , 

T°(e; e p ( 4 _i)) = { (4.15) if q tr ( 0) > 0 , q tr (V6f c g tr /2Gm 0 ) < 0 , 

O if q tr (\/6f c g tr /2Gm 0 ) > 0. 


(4.16) 


5 An abstract model 


The aim of this section is an extension of Theorem |3.1 and 4T on a specific class of elastoplastic 
models that are usually formulated in the Haigh-Westergaard coordinates. We consider an abstract 
model containing the isotropic hardening and the plastic flow pseudo-potential. Given the history of 
the strain tensor s = e(t), t e [t 0 , t max ], and the initial values 

s p (t 0 ) = s p , s p (t 0 ) = e p 0 . 

Find the generalized stress (cr(t), n(t)) and the generalized strain (s p (t),£ p (t)) such that 


s = s e + s p , 

cr = O e : s e , k — H(£ p ), 

s p G Xd a g(a,n), > 

if? = X£(cr, hi), 

A > 0, /(cr, hi) < 0, A/(cr, n) — 0. , 


hold for each instant t e [4 Umax] ■ 

Further, we have the following assumptions on ingredients of the model: 

1 . f(er,n) = f(p(cr), g(cr), g e (cr), k), where / is increasing with respect to g and g, convex and 
continuously differentiable at least in vicinity of the yield surface. 

2 . g{cr,hi) = g(p(cr), g(er), n), where g is an increasing function with respect to g, convex and 
twice continuously differentiable at least in vicinity of the yield surface. 
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3. H is a nondecreasing, continuous and strongly semismooth function satisfying H( 0) = 0. 

4. £(cr , k) = £(p(cr), g(cr), g(cr), n) is a positive function. 

5. Invariants p, g, g e and g are the same as in Section [2] 

Notice that the assumptions on / and g guarantee convexity of / and g using properties of r e 
introduced in [39]. Let gv(p, g) '■= 9g(Pjg) Then one can write the plastic flow rule as 

follows: 

e p = X [g v (p, g)I /3 + g e {jp, e)n\, n e dg(cr). 

The fc-th step of the incremental constitutive problem received by the implicit Euler method reads 
as follows. Given e := e(tk), <J tr O e : (e(Lt) — £ p (tk- 1 )) and := eFiftk- 1 ). Find <j = cr{tk), 
s? = tk) and AX satisfying: 

cr = cr tr - AA [Kg v (p(cr), g{cr))I + 2Gg e (p(cr), g{cr))n] , n G 9p(cr), 1 

£* = £*>.*»• + AA%,k), > (5.1) 

AA > 0 , f(tr,H{ee))< 0, AA/(<r, = 0. I 


If we use the elastic predictor/plastic corrector method then we derive the following straightfor¬ 
ward extension of Theorem |3.1| and Theorem |4.1| within the plastic correction. 


Theorem 5.1. Let f(cr tr ,H(e p ' tr )) >0. If(cr,e p , AA) is a solution to problem (5.1) then (p, g,^, AA), 
AA > 0 , is a solution to the following system: 


P = P tr ~ AA Kg v (p, g), 
g = [g tr - AA 2 Gg e (p, p)] + , 
eP = £P’ tr + A A l (p, g, gr( cos 6 tr )), 
/ (p, g, gr e (cosd tr ), H(eP)) = 0 . 


Conversely, if (p, g, k , AA) is a solution to (5.2) then (<r, k, AA) solves (5.1), where 


cr = 


pi + pn* r if g > 0 , 

pi if 6 = 0 . 


(5.2) 


(5.3) 


Notice that it is generally impossible to a priori decide about the type of the return as in the 
models introduced above. To be in accordance with the current approach introduced e.g. in na one 


can split (5.2) into the following two systems: 


p + AXKg v (p,0) 

= p tr 

£P- AAl(p, 0 , 0 ) 

= £P’ tr 

/ (p, 0,0, Hie 9 )) 

= 0 

p + AXKg v (p, g) 

tv 

= p 

g + AX2Gg e (p, g) 

= g tr 

?? — AX £ (p, g, gr( cos 6 tr )) 

= £P’ tr 

f(p, 6, gr e (cos0 tr ),H(e p )) 

= 0 


for g = 0 ( return to the apex (apices)), 


(5.4) 


for q > 0 ( return to the smooth portion), (5.5) 


and guess which of these systems provides an admissible solution. Beside the blind guessing, the 


current approach has another drawback: it can happen that (5.2) has a unique solution and mutually 


one of the systems (5.4), (5.5) does not have any solution or have more than one solution. Therefore, 


we recommend to solve (|5.2|) directly by a nonsmooth version of the Newton method with the standard 
initial choice p° = p tr , 


g° = g tr , 


K 


0 = n tr and AA° = 0 . 
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6 Numerical realization of the incremental boundary value 
elastoplastic problem 

Consider an elasto-plastic body occupying a bounded domain QCK 3 with the Lipschitz continuous 
boundary T. It is assumed that T = T d U Tjv, where T o and T n are open and disjoint sets. On 
T £>, the homogeneous Dirichlet boundary condition is prescribed. Surface tractions of density f t are 
applied on T n and the body is subject to a volume force f v . 

Notice that the above defined stress, strain and hardening variables depend on the spatial vari¬ 
able x E Q, i.e. cr k = cr k (x), etc. Let V := {v e [iL 1 (0)] 3 | v = 0 on T#} denote the space of 
kinematically admissible displacements. Under the infinitesimal small strain assumption, we have 
e(v) = \ (Vu + (Vu) T ) , v e V. 

Substitution of the stress-strain operator T into the principle of the virtual work leads to the 
following problem at the k-th step: 

(Pk) find u k e V : j T (e(w fc ); e£_j, e^_j) : e(v) dft = [ f Vk .vdFt + f f t>k .v dT \/v € V, 

J Yl J Yl J Y n 

where f Vk and f tk are the prescribed volume, and surface forces at t k , respectively. After finding 
a solution u k , the remaining unknown fields £ k ,£ k important for the next step can be computed at 
the level of integration points. Problem (P k ) can be standardly written as the operator equation in 
the dual space V to V: J r k(uk) = £ k , where 

(■ P k (u),v) = [ T(e(i4);e^_ 1 ,^_ 1 ) : e(v) dfi V«,r6V, 

Jn 

(h,v) = I f Vtk .vdU+ [ f tjk .vdT Wv e V. 

J Yl J Tjv 

Since we plan to use the semismooth Newton method, we also introduce the operator K, k : V —* 
£(V,V / ) as follows: 

(IC k (u)v,w) — I T° e(v) : e(w) dQ Vu,v,weV. 

Jn 

To discretize the problem in space we use the finite element method. Then the space V is 
approximated by a finite dimensional one, V/,,. If linear simplicial elements are not used then it 
is also necessary to consider a suitable numerical quadrature on each element. Let Pk,h, K-k,h, tk,h 
denote the approximation of operators P k , JC k , kk, respectively, and F k : M n —* M n , K k : W 1 —> M nxn , 
l k E M n be their algebraic counterparts. Then the discretization of problem (P k ) leads to the system 
of nonlinear equations, F k (u k ) = l k , and the semismooth Newton method reads as follows: 
Algorithm 1 (Semismooth Newton method). 

1 : initialization: u k = u k _\ 

2 : for i — 0 , 1 , 2 ,.. . do 

3: find 5u' 1 G V : K k (ul)Su l = l k — F k (u\) 

4: compute u t k 1 = u\. + hu 1 

5: if ||5w i ||/(||uj, +1 || + ||wj.||) < rNewton then stop 

6: end for 
7: set u k = u j, +1 . 

If T is strongly semismooth in then F k is strongly semismooth in M n . Notice that the 

strong semismoothness is an essential assumption for local quadratic convergence of this algorithm. 
In numerical examples introduced below, we observe local quadratic convergence when the tolerance 
is sufficiently small. In particular, we set eNewtcm = 10 -12 . 
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7 Numerical example - slope stability 

The improved return-mapping schemes in combination with the semismooth Newton method have 
been partially implemented in codes SIFEL [Tj and MatSol [21] • Here, for the sake of simplicity, we 
consider the slope stability benchmark HD Page 351] for the presented models, the Drucker-Prager 
(DP) and the Jirasek-Grassl (JG) ones. The benchmark is formulated as a plane strain problem. 
We focus on: a) incremental limit analysis and b) dependence of loading paths on element types and 
mesh density. For purposes of such an experiment, special MatLab codes have been prepared to be 
transparent. These experimental codes are available in [2] together with selected graphical outputs. 

The geometry of the body is depicted in Figure [2] or [3] The slope height is 10 m and its inclination 
is 45°. On the bottom, we assume that the body is fixed and, on the left and right sides, zero normal 
displacements are prescribed. The body is subjected to self-weight. We set the specific weight 
pg = 20kN/m 3 with p being the mass density and g the gravitational acceleration. Such a volume 
force is multiplied by a scalar factor, (. The loading process starts from ( = 0. The gravity load 
factor, (, is then increased gradually until collapse occurs. The initial increment of the factor is set 
to 0.1. To illustrate loading responses we compute settlement at the corner point A on the top of 
the slope depending on (. 

As in [TTJ Page 351], we set E = 20 000kPa, u = 0.49, (ft = 20° and c = 50kPa, where c denotes 
the cohesion for the perfect plastic model. Hence, G = 67114kPa and K = 3 333 333 kPa. In 
comparison to HU, we use the presented models instead of the Mohr-Coulomb one. The remaining 
parameters for these models will be introduced below. 

We analyze the problem for linear triangular (PI) elements and eight-pointed quadrilateral (Q2) 
elements. In the latter case, (3 x 3)-point Gauss quadrature is used. For each element type, a 
hierarchy of four meshes with different densities is considered. The Pl-meshes contain 3210, 12298, 
48126, and 190121 nodal points, respectively. The Q2-meshes consist of 627, 2405, 9417, and 37265 
nodal points, respectively. The coarsest meshes for PI and Q2 elements are depicted in Figure [2] 
and[3j Let us complete that the mesh in Figure [2] is uniform in vicinity of the slope and consists of 
right isoscales triangles with the same diagonal orientation. Further, it is worth mentioning that the 
Pl-meshes are chosen much more finer in vicinity of the slope than their (^2-counterparts within the 
same level. 



Figure 2: The coarsest mesh for PI elements. 


Figure 3: The coarsest mesh for Q2 elements. 


7.1 The Drucker-Prager model 

The Drucker-Prager parameters r/, rj and £ are computed from the friction angle, (ft , and the dilatancy 
angle, ift, as follows CU: 

3tan0 _ 3tan-0 3 

\J 9 T 12(tan (ft) 2 ’ ^ y/9 + 12(tan£>) 2 ’ y/9 + 12(tan (ft) 2 
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At first, we introduce results obtained for the model with associative perfect plasticity. In such a case, 
tj} = (p, Co = c and H — 0 kPa. The received loading curves for the investigated meshes and elements 
are depicted in Figure [| and [5j Although Pl-meshes are much finer, we observe more significant 
dependence of the curves on the mesh density for P 1-elements than for (32-elements. Also computed 
limit load factors are greater and tend more slowly to a certain value as h —y 0 + for Pl-meshes than 
for (32-meshes. The expected limit value is 4.045 as follows from considerations introduced in [7j. 
Using the finest PI and Q2 meshes, we receive the values 4.230, and 4.056, respectively. 

In general, higher order elements are recommended when a locking effect is expected. In this 
example, it can be caused due to the presence of the limit load and/or v « 1/2. On the other hand, 
the strong dependence on mesh density is influenced by other factors like mesh structure or choice of 
a model. For example, this dependence is not so significant for the Jirasek-Grassl model (see the next 
subsection). Further, in [T9], there is theoretically justified and illustrated that the dependence of 
the limit load on the mesh density is minimal for bounded yield surfaces and that an approximation 
of unbounded yield surfaces by bounded ones (the truncation) leads to a lower bound of the limit 
load. 




Figure 4: Loading paths for the associative perfect Figure 5: Loading paths for the associative perfect 

plastic model and PI elements. plastic model and Q 2 elements. 

For illustration, we add Figure [6] and [7] with plastic multipliers and total displacements at collapse, 
respectively. The figures are in accordance with literature. 



Figure 6: Plastic multipliers at collapse for the Figure 7: Displacements at collapse (detail) for 
finest <32-mesh. the (32-mesh with 2405 nodes. 

To compare the current return-mapping scheme with improved one, we have also considered the 
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nonassociative model with nonlinear hardening where 


ip = 10°, c 0 = 40kPa, H = 10000 kPa, //(A 5 ) 


min 


c — c 0 , 


He?- 


H 2 

4(c - c 0 ) 



Here, H represents the initial slope of the hardening function and the material response is perfect 
plastic for sufficiently large values of the hardening variable. This nonassociative model yields a 
slightly lower values of the limit load factors and also the other results are very similar to the 
associative model. The related graphical outputs are available in [21 SS-DP-NH], Further, in vicinity 
of the limit load, we have observed lower rounding errors for the improved return-mapping scheme 
and thus lower number of Newton steps is necessary to receive the prescribe tolerance than for the 
current scheme. However, the computational time for both schemes are practically the same since 
return to the apex happens only on a few elements lying in vicinity of the yield surface. 


7.2 The simplified Jirasek-Grassl model 

To be the simplified Jirasek-Grassl (JG) model applicable for the investigated soil material we fit 
its parameters using the associative perfect plastic Drucker-Prager (DP) model as follows: e = 1, 
f c = 3c^/(v / 3 — ?/), ft = 0 , B g = 1000, s = 5, A g = sr] and m 0 = — 6. Recall that e = 1 implies 

g e = g. Further the value of f c corresponds to the uniaxial compressive strength computed from 
the Drucker-Prager model. To eliminate the influence of the exponential term in the function m g) 
the value of B g is chosen sufficiently large. Then the model is insensitive on f t and one can vanish 
it. Finally, we require the same flow direction for both the models under the uniaxial compressive 
strength. Since the yield function in the JG model is normalized in comparison to the DP model it 
is convenient to introduce the following relation between the plastic multipliers: AA dp = ■f-AAjc', 

JG 

where s is a scale factor. Then the values of m 0 and A g are determined from the following equations: 

SV = m' g (f c / 3 ) « A g , -)= = VQ + ^=. 


To be mo positive, s must be greater than 2x/3. To be in accordance with results of the DP model, 
we set s = 4. 90 Comparison of yield surfaces (in the meridean plane) and flow directions for the 
DP and JG models is illustrated in Figure [ 8 } Here, the fixed value AA dp = 0.001 is used for vectors 
representing the flow directions. 

Loading curves for the investigated PI, Q 2 meshes and the JG model are depicted in Figure[9]and 
We observe much faster convergence of the Pl-loading curves than for the DP model. Moreover, 
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the results for PI and Q2 elements are comparable. The computed values of the limit load factor on 
the finest PI and Q2 meshes are 4.124, and 4.107, respectively. 


8 Conclusion 

The main idea of this paper is that the sub differential formulation of the plastic flow rule is also 
useful for computational purposes and numerical analysis. Namely, it has been shown that such an 
approach improves the implicit return-mapping scheme for non-smooth plastic pseudo-potentials as 
follows. 

• The unique system of nonlinear equations is solved regardless on a type of the return. 

1 For smaller values of s, the limit load factor is underestimated and for greater values of s the limit load factor is 
overestimated. 
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Figure 8: Comparison of the yield surfaces (in the meridean plane) and the flow directions for the DP 
model (black) and the JG model (grey). 




Figure 9: PI - loading paths for the simplified Figure 10: Q2 - loading paths for the simplified 
Jirasek-Grassl model. Jirasek-Grassl model. 
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• It can be a priori determined the type of the return from a given trial state for some models 
(without knowledge of the solution). 

• The scheme can be more correct than the current one, and its form enables to study properties 
of constitutive operators like existence, uniqueness and semismoothness. 

In this paper (PART I), the new technique has been systematically built on a specific class of models 
containing singularities only along the hydrostatic axis. Beside an abstract model, two particular 
models have been studied: The Drucker-Prager and the simplified Jirasek-Grassl model. However, 
the presented idea seems to be more universal. For example, it has been successfully used for the 
Mohr-Coulomb model in ’’PART II” [57] . 
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